
packages = c("brglm","ggcorrplot","interplot", "stargazer", "tidyverse",
             "here","haven","mice","modelsummary", "sandwich", "lmtest",
             "MASS", "ROCR", "generalhoslem", "brant", "bcROCsurface", "gridExtra",
             "erer", "margins", "interplot.medline")
#Warning! This will install R packages in your system if not yet installed. 
#load or install & load all packages 
package.check <- lapply(
  packages,
  FUN = function(x) {
    if (!require(x, character.only = TRUE)) {
      install.packages(x, dependencies = TRUE)
      library(x, character.only = TRUE)
    }
  }
)


GH_SS_Replication <- read_dta(here::here("GH_SS_Replication.dta"))
VDem <- read.csv(here::here("Vdem","V-Dem-CY-Full+Others-v12.csv"))

 GH_SS_clean <- GH_SS_Replication %>% 
  subset(select = c(warnumber, warname, battlename, year, ccode, participate,
                    win, aggress, msadopt,sidebmscat2,msadoptcompare, troopsengaged, opptroopsengaged,
                    terrain, strat, troopbalance, sideAbof, cinc, 
                    bankseducation,bankseduccompare,bankscoups, bankscoupscompare, polity, politycompare,
                    politysquared, politysideb, politysidebsquared)) %>%
  dplyr::rename(COWcode = ccode) %>% 
  mutate(deffender = ifelse(aggress ==0,1,0),
         pre_war_year = year -1) %>% drop_na(warnumber, battlename)



vdem_clean <- VDem  %>% subset(select = c(country_name,
                                          COWcode,
                                          year,
                                          v2x_corr,
                                          e_gdppc,
                                          v2x_ex_military,
                                          v2x_libdem,
                                          v2x_polyarchy,
                                          v2regsupgroups_6)) 

vdem_clean$pre_war_year <-vdem_clean$year -1 
vdem_clean$COWcode <- ifelse(vdem_clean$country_name == "Austria", 300, vdem_clean$COWcode )
vdem_clean$COWcode <- ifelse(vdem_clean$country_name == "Piedmont-Sardinia", 325, vdem_clean$COWcode )


growth_rate <- function(x)(x/lag(x)-1)*100 
vdem_clean <- vdem_clean %>% group_by(COWcode) %>% mutate(gdppc_growth = growth_rate(e_gdppc)) 
GH_SS_data <- left_join(GH_SS_clean,vdem_clean, by = c("COWcode", "pre_war_year")) %>% subset(select = c(-year.y)) %>% dplyr::rename(year = year.x)


#######create rel corr ########
GH_SS_data$v2x_corr_init <-GH_SS_data$v2x_corr*GH_SS_data$aggress
GH_SS_data$v2x_corr_target <- GH_SS_data$v2x_corr*GH_SS_data$deffender

GH_SS_data3 <- GH_SS_data %>% dplyr::select(warnumber,
                                            battlename,
                                            aggress,
                                            deffender,
                                            v2x_corr_init,
                                            v2x_corr_target) 

GH_SS_data2 <-GH_SS_data %>% group_by(warnumber,battlename, aggress) %>% summarise(init_corr_mean = mean(v2x_corr_init, na.rm = TRUE),
                                                                        target_corr_mean = mean(v2x_corr_target, na.rm = TRUE))

GH_SS_data2 <-GH_SS_data2 %>% group_by(warnumber, battlename) %>% summarise(init_corr_mean = sum(init_corr_mean),
                                                                    target_corr_mean = sum(target_corr_mean) )

GH_SS_data<- left_join(GH_SS_data, GH_SS_data2, by = c("warnumber", "battlename"))
GH_SS_data$opp_corr <- ifelse(GH_SS_data$aggress==1, GH_SS_data$target_corr_mean, GH_SS_data$init_corr_mean)
GH_SS_data$relative_corr <- GH_SS_data$v2x_corr - GH_SS_data$opp_corr

#Check Percentage of missing data for key variables 
GH_SS_data <-GH_SS_data  %>% filter(participate != 0)
test <-GH_SS_data  %>% filter(is.na(participate))

pMiss <- function(x){sum(is.na(x))/length(x)*100}
apply(GH_SS_clean,2,pMiss)
apply(GH_SS_data,2,pMiss)
apply(GH_SS_Replication,2,pMiss)


#Models from Replication of Grauer and Horowitz (2012)
#*/ M4 */
#  probit win msadopt polity politysquared cinc strat terrain troopsengaged opptroopsengaged bankseducation loggdppercapita bankscoups, cluster(ccode)

#*/ M5 */
#  probit win msadoptcompare polity politysquared politysideb politysidebsquared sideAbof terrain strat troopbalance bankseduccompare bankscoupscompare, cluster(ccode)


m1 <- glm(win ~ relative_corr + msadopt,
         data = GH_SS_data,  family = binomial(link = "probit"))
summary(m1)
m1sr<- coeftest(m1, 
                vcov = sandwich(x =m1, adjust = TRUE),
                cluster = COWcode) #stata robust, slightly different because of of PC rather than HC

m1a <- glm(win ~ relative_corr + msadopt + polity + politysquared + cinc+ strat + terrain + troopsengaged + opptroopsengaged + e_gdppc + bankseducation + bankscoups ,
           data = GH_SS_data,  family = binomial(link = "probit"))
summary(m1a)
m1asr<- coeftest(m1a, 
                 vcov = sandwich(x =m1a, adjust = TRUE),
                 cluster = COWcode) #stata robust, slightly different because of of PC rather than HC

m1b <- glm(win ~ relative_corr + msadopt + v2x_libdem + troopsengaged + cinc + strat + terrain,
          data = GH_SS_data,  family = binomial(link = "probit"))
summary(m1b)
m1asr<- coeftest(m1b, 
                vcov = sandwich(x =m1b, adjust = TRUE),
                cluster = COWcode) #stata robust, slightly different because of of PC rather than HC

stargazer(m1, m1a, m1b, 
          se = list(m1asr[,2], m1asr[,2],m1asr[,2]),
          covariate.labels = c("Relative Corruption", "Modern System Adoption Level", "Polity Score", 
                               "Polity Squared","Liberal Democracy", "Material Capabilities",
                               "Strategy", "Terrain", "Troops Engaged",
                               "Opponent Troops Engaged", "GDPPC","Human Capital", "Civil-Military Relations"
                               ))

m2 <- glm(win ~ relative_corr + msadoptcompare,
          data = GH_SS_data,  family = binomial(link = "probit"))
summary(m2)
m2sr<- coeftest(m2, 
                vcov = sandwich(x =m2, adjust = TRUE),
                cluster = COWcode) #stata robust, slightly different because of of PC rather than HC

m2a <- glm(win ~ relative_corr + msadoptcompare + polity + politysquared +  politysideb + politysidebsquared + 
               sideAbof + strat + terrain + troopbalance + bankseduccompare +  bankscoupscompare,
             data = GH_SS_data,  family = binomial(link = "probit"))
summary(m2a)
m2asr<- coeftest(m2a, 
                 vcov = sandwich(x =m2a, adjust = TRUE),
                 cluster = COWcode) #stata robust, slightly different because of of PC rather than HC


m2b <- glm(win ~ relative_corr + msadoptcompare + v2x_libdem + sideAbof +troopbalance + strat + terrain,
           data = GH_SS_data,  family = binomial(link = "probit"))
summary(m2b)
m2bsr<- coeftest(m2b, 
                 vcov = sandwich(x =m2b, adjust = TRUE),
                 cluster = COWcode) #stata robust, slightly different because of of PC rather than HC
 

stargazer(m2, m2a, m2b, 
          se = list(m2sr[,2], m2asr[,2],m2bsr[,2]),
          covariate.labels = c("Relative Corruption", "Comparative Modern System Adoption Level", "Polity Score", 
                               "Polity Squared","Polity Side B Score", "Polity Side B Squared ", "Liberal Democracy", 
                               "Comparative Material Capabilities", "Strategy", "Terrain", "Comparative Troops Engaged",
                               "Comparative Human Capital", "Comparative Civil-Military Relations"
          ))

m3 <- lm(msadopt ~ v2x_corr ,
           data = GH_SS_data,  family = binomial(link = "probit"))
summary(m3)
m3sr<- coeftest(m3, 
                 vcov = sandwich(x =m3, adjust = TRUE),
                 cluster = war) #stata robust, slightly different because of of PC rather than HC
 
 
m3a <- lm(msadopt ~ v2x_corr  + v2x_libdem +bankseducation + bankscoups + cinc +  e_gdppc + year,
            data = GH_SS_data,  family = binomial(link = "probit"))
summary(m3a)
m3asr<- coeftest(m3a, 
                vcov = sandwich(x =m3a, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC


GH_SS_data$binary_msadopt <- ifelse(GH_SS_data$msadopt > 1,1,0)

m4 <- glm(binary_msadopt ~ v2x_corr,
           data = GH_SS_data,  family = binomial(link = "probit"))
summary(m4)
m4sr<- coeftest(m4, 
                vcov = sandwich(x =m4, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC


m4a <- glm(binary_msadopt ~v2x_corr  + v2x_libdem + bankseducation + bankscoups + cinc +  e_gdppc + year,
           data = GH_SS_data,  family = binomial(link = "probit"))
summary(m4a)
m4asr<- coeftest(m4a, 
                vcov = sandwich(x =m4a, adjust = TRUE),
                cluster = war) #stata robust, slightly different because of of PC rather than HC

stargazer(m3, m3a, m4, m4a, 
          se = list(m3sr[,2], m3asr[,2],m4sr[,2],m4asr[,2]),
          covariate.labels = c("Relative Corruption", "Liberal Democracy", 
                              "Human Capital", " Civil-Military Relations", "Material Capabilities", "GDPPC", "Year"),
          dep.var.labels= c("Modern System Adoption Level ",
                            "Binary Modern System Adoption")
          )


